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3D simulations of wet foam coarsening evidence a self similar growth regime 
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Abstract 

In wet liquid foams, slow diffusion of gas through bubble walls changes bubble pressure, volume and wall curvature. Large 
bubbles grow at the expenses of smaller ones. The smaller the bubble, the faster it shrinks. As the number of bubbles in a given 
volume decreases in time, the average bubble size increases: i.e. the foam coarsens. During coarsening, bubbles also move relative 
to each other, changing bubble topology and shape, while liquid moves within the regions separating the bubbles. Analyzing 
the combined effects of these mechanisms requires examining a volume with enough bubbles to provide appropriate statistics 
throughout coarsening. Using a Cellular Potts model, we simulate these mechanisms during the evolution of three-dimensional 
foams with wetnesses of cp = 0.00, 0.05 and 0.20. We represent the liquid phase as an ensemble of many small fluid particles, 
which allows us to monitor liquid flow in the region between bubbles. The simulations begin with 2x10^ bubbles for 0 = 0.00 and 
1.25 X 10^ bubbles for 0 = 0.05 and 0.20, allowing us to track the distribution functions for bubble size, topology and growth rate 
over two and a half decades of volume change. All simulations eventually reach a self-similar growth regime, with the distribution 
functions time independent and the number of bubbles decreasing with time as a power law whose exponent depends on the wetness. 


1. Introduction 

Liquid foams include soap froths and food and industrial 
foams, as well as many liquid-phase pollutants. Foams consist 
of gas bubbles surrounded by a continuous liquid phase. Bub¬ 
ble walls are very thin; their thickness strongly depends on the 
intensity of mutual repulsion between surfactant-loaded gas- 
liquid interfaces and only weakly on foam wetness (0=liquid 
volume/total volume) below a critical wetness at which the 
foam transforms into a bubbly liquid. When three bubble walls 
meet, they form edges, known as Plateau borders, which are 
tube-like structures with a cross section which is roughly a tri¬ 
angle with its straight sides replaced by concave circular arcs. 
Wet foams and dry foams both have bubbles separated by thin 
walls, the difference being that wetter foams have larger Plateau 
borders mm. 
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Slow diffusion of gas through bubble walls and Plateau bor¬ 
ders, due to the Laplace pressure difference, changes bubble 
pressure, volume and wall curvature. Bubbles, Plateau bor¬ 
ders and walls move and relax towards mechanical equilibrium, 
while liquid flows inside the Plateau borders and walls as well. 
Bubbles larger than the mean bubble size (or bubbles with more 
than the mean number of faces for dry foams) grow at the ex¬ 
penses of smaller (fewer-faced) ones. Smaller (fewer-faced) 
bubbles shrink more and more quickly, and eventually disap¬ 
pear. The number of bubbles decreases in time, and the average 
bubble size increases: the foam coarsens dElEl (figure [^. 

X-ray tomography experiments ID and numerical simula¬ 
tions nni have found that dry foams reach a self-similar 
growth regime where the topology and relative bubble size dis¬ 
tributions are steady while the mean bubble radius increases 
as as expected from dimensional analysis (HI. Hence, the 
mean bubble volume increases as and the number of bub- 
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Figure 1: Snapshots of a 3D foam-coarsening experiment. First (left) and last 
(right) images of a 3D movie obtained by (4), in a foam with 0 = 0.17. 

bles decreases as In the opposite limit, in a bubbly liq¬ 

uid, gas flows from small bubbles to the liquid phase, where 
it diffuses, to eventually flow back into larger bubbles due to 
Laplace pressure difference, now between the bubbles’ inter¬ 
nal pressures and the gas saturation pressure in the liquid. This 
Ostwald ripening dynamics also generates a self-similar growth 
regime, where the relative bubble size distribution is constant, 
and the mean bubble radius grows as 0 . Hence, the mean 
bubble volume grows as t and the number of bubbles decreases 
as 

To investigate coarsening for foams of intermediate wet¬ 
nesses, Isert, Maret and Aegerter Col levitated an aqueous 
foam in a magnetic held: the water’s diamagnetism counterbal¬ 
anced gravity, so the wetness remained uniform throughout the 
foam. At long times, for foams with wetness less than 0 = 0.25 
the mean bubble radius grew as as in a dry foam, for foams 
with wetness greater than 0 = 0.35 the mean bubble radius 
grew as as in a bubbly liquid, and the mean bubble radius 
growth exponent decreased from 1/2 to 1/3 between 0 = 0.25 
and 0 = 0.35. Determining whether power law growth in bub¬ 
ble size at these wetnesses reflected self-similar growth with 
constant relative bubble size distributions would require deter¬ 
mining individual bubble properties within the levitated foam. 
Such a measurement would be quite difficult, though it might 
be possible using e.g. X-ray tomography fdj. 

Numerical simulations provide access to properties of in¬ 



Figure 2: Representation of the grid neighborhood for Potts energy calculations, 
up to the 5^^ layer, summing over 56 neighbors. 

dividual bubbles (volume, face-number and edge-number dis¬ 
tributions, and their spatial correlations) with arbitrary spatial 
and temporal resolution. Complicating such simulations are the 
multiple mechanisms, length and time scales affecting coarsen¬ 
ing. Slow volume changes in bubbles cause the gradual increase 
of Plateau border cross section. Both bubble and Plateau-border 
growth can induce topological transformations in bubble and 
Plateau-border adjacency, such as neighbors swapping (Tl), 
bubble disappearance (T2), or Plateau-border fusion, which can 
cause further rapid cascades of rearrangement. These rear¬ 
rangements, in turn, change the slower dynamics of gas ex¬ 
change (between neighboring bubbles, and between bubbles 
and the liquid phase) and flow in Plateau borders. Length 
scales include the thickness of the thin walls (which can range 
from hundreds of nanometers to many microns), the size of the 
Plateau borders (which can range from hundreds of microns to 
millimeters), the lengths of Plateau borders and the radii of bub¬ 
bles (which can range from sub-millimeter to centimeters) de¬ 
pending on the foam, wetness and degree of coarsening. The 
close feedback among all these mechanisms and scales makes 
simulation of wet foam dynamics challenging ifTTIl . 

The Cellular Potts model (CPM) can simulate the structure 
of foams in mechanical equilibrium, with a resolution (number 
of pixels per bubble) which matches or exceeds that of experi¬ 
mental images in two or three dimensions. The CPM can also 
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simulate the quasi-static dynamics of foam coarsening, as slow 
gas diffusion leads to a succession of mechanically equilibrated 
foam structures. The CPM correctly reproduces bubble rear¬ 
rangements and disappearance. 

The CPM developed from Potts model simulations of coars¬ 
ening in metallic polycrystals ca, because the mass flux 
between bubbles and crystals both obey Pick’s law (a local 
dynamics which produces the same growth dynamics as the 
Young-Laplace law for bubbles). Such CPM simulations agree 
in almost all respects with experimental observations of dry liq¬ 
uid foams na. However, in the CPM, mass flow between bub¬ 
bles and surface-tension-driven relaxation of wall shape occur 
over similar time scales, while the former is much slower than 
the latter in almost all liquid foams. Extending the CPM method 
by interleaving steps in which mass is free to diffuse between 
bubbles with steps in which we constrain the bubble volumes to 
their current values but allow the bubble walls to relax allows 
us to reproduce the experimental hierarchy of time scales. The 
more relaxation steps per diffusion step, the faster the relaxation 
rate relative to the diffusion rate 0. 

CPM simulations of 3D dry foam coarsening (61 [71 and 
2D wet foam coarsening both reach self-similar growth 
regimes. Here we compare 3D CPM simulations (with ad¬ 
ditional shape relaxation steps) of foams with wetnesses of 
0 = 0.00, 0.05 and 0.20. All three reach self-similar growth 
regimes with enough bubbles to yield statistically significant 
distribution functions. 

2. Methods 

2.7. Foam structure 

Like experimental images, CPM 3D foam simulations rep¬ 
resent space as a cubic grid of voxels, and each bubble as a 
connected set of voxels with particular properties. We employ 
a space of Vt = 200^ voxels with periodic boundary conditions. 



Figure 3: Cross-sections through a 200^ voxels 3D CPM simulation of a liquid 
foam with periodic boundary conditions, for 0 ^ 0.05 after 5600 MCS. Left: 
without interface relaxation. Right: with interface relaxation. 

initially filled with 2 x 10^ bubbles for 0 = 0.00 and 1.25 x 10^ 
bubbles for 0 = 0.05 and 0.20. Each voxel, located at a posi¬ 
tion r, has a label S , corresponding to a bubble or portion of 
liquid. Defining the liquid phase as a single domain would al¬ 
low liquid to move instantaneously within the foam, effectively 
giving it zero inertia and viscosity. To ensure local conserva¬ 
tion of liquid volume and to keep the liquid velocity finite, we 
follow Eortuna et al C3, and subdivide the liquid phase into 
numerous small fluid particle subdomains (see figure S1 in the 
supplementary materials). 

The CPM energy function includes a surface energy and a 
volume constraint: 

£ = 7(5^,+ L ^(>5) (Vs - (1) 

r v{r) S 

where the sum over v(F} extends through the 5^^ neigh¬ 
bors around the voxel at r (56 neighbors), to reduce lattice- 
anisotropy effects ca, and is only evaluated for neighboring 
voxels with distinct labels S . Eigure [^illustrates the grid neigh¬ 
borhood in the energy calculation. The second term in the r.h.s. 
of Eq. 0 constrains the volume of the labeled region S, Vs, 
around its target volume, when /1(5'), which is the com¬ 

pression modulus (inverse compressibility) of the material, is 
> 0 . 

While an “ideal” fluid particle would be negligibly small 
and perfectly incompressible, we assign our fluid particles a 
ytarget _ y yoxels SO they are much smaller than typical bub- 
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Figure 4: Images of a CPM simulation of a liquid foam with 0 = 0.20, run in a 
grid with 300^ voxels. Initial (left) and final (right) states. 

bles, but large enough to move freely (see figure S2 in the sup¬ 
plementary materials for the effect of varying the target volume 
of the fluid particles, which affects coarsening dynamics). The 
compression modulus Xuq = 10 is not critical, since coarsening 
dynamics are robust under moderate variation of this parameter 
(see figure S3 in the supplementary materials). However, mak¬ 
ing Aiiq too large {i.e. the compressibility too small) leads to 
non-physical “freezing” of the liquid in the Plateau borders. 

Simulations begin with a number of fluid particles randomly 
distributed over the grid. The foam wetness 0 is simply the 
number of voxels in fluid particles (controlled by their number 
and target volume) divided by the lattice size. We All the rest 
of the grid with bubbles with initial target volumes randomly 
drawn from a log-normal, normal or bidisperse probability dis¬ 
tribution. 

The first term of the r.h.s. of Eq. Q defines the interfacial 
energy between the different types of regions: 

JQiqJiq) = 0.10, 

J(liq, gas) = J(gas,liq) = 1.00, 

J{gas, gas) = 1.99. (2) 

J{liq, liq) would be zero in a real liquid, but a small positive 
value helps to maintain the fluid particles as connected domains 
without significantly affecting their dynamics. J{liq, gas) is the 
energy per unit area of the gas-liquid interface and J(gas, gas) 


CD 


(D 




Figure 5: Growth rate. Upper panel: Growth rate as function of the normalized 
bubble volume V/ (V) for simulations starting from a log-normal bubble target- 
volume distribution with 0 = 0.20 (top). The black line shows the average over 
10 simulation runs for 31 time points after the onset of self-similar growth. 
Lower panel: Self-similar growth rate as a function of the normalized bubble 
volume V/ {V) for different wetnesses. While similar, these growth rates are not 
identical, showing a significant dependence on 0. 


is of the order of 2J{liq, gas) because a wall separating two bub¬ 
bles consists of two gas-liquid interfaces. If J(gas, gas) were 
strictly equal to 2J(liq, gas), the fluid particles would spread in 
the walls as in a bubbly liquid. We set J(gas, gas) slightly lower 
than 2J(liq, gas) to account for the positive disjoining pressure 
{2J{liq, gas) - J{gas, gas)) between the soap Aims in the walls, 
which serves to keep the walls thin, as observed in experiment. 
The greater the disjoining pressure, the greater the fraction of 
liquid that localizes in the Plateau borders and the less that pen¬ 
etrates into the bubbles walls. As a result, we expect coarsen¬ 
ing to be faster for higher disjoining pressures. For the range 
2J(liq, gas) - J(gas, gas) = 0.1, 0.01 and 0.001, the simulation 
results are robust and the liquid’s spatial distribution is essen¬ 
tially the same (see flgure S4 in the supplementary materials). 
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The wall thickness and Plateau-border shape depend on the 
interaction range (v(f) in Eq. (1)) and the size of the fluid parti¬ 
cles. 

2.2. Foam dynamics 

Conflgurations change stochastically in the CPM. We ran¬ 
domly pick a pair of first neighbor voxels with different la¬ 
bels and calculate the current energy, Ei, and the energy after 
copying the label S from the first to the second voxel, Ef. If 
AE” = Ef - Ef < 0 we always make the copy. If AE > 0 
we make the copy with probability exp(-A£'/r), where T is 
a Boltzmann-like fluctuation amplitude which we choose just 
large enough to overcome lattice pinning effects (the macro¬ 
scopic fluctuation amplitude in a real foam would be close to 
0, but real foams do not experience lattice pinning). The simu¬ 
lation unit which corresponds to time is the Monte Carlo Step 
(MCS): 1 MCS consists of N voxel-copy attempts, where N is 
the total number of voxels in the grid. 

Like all models, the CPM’s simplification of reality may in¬ 
troduce artifacts. The CPM models mass transfer as voxel trans¬ 
fer. The transfer of mass between two bubbles in contact across 
a thin wall obeys Pick’s law in both CPM and experiment, so 
the key question is the rate of transfer as a function of parame¬ 
ters. In both experiment and simulation, we expect the transfer 
of mass across the Plateau borders to be slower than across thin 
walls. However, the dependence of the rate of transfer of bub¬ 
ble mass between bubbles via the Plateau borders depends less 
obviously both on model parameters and geometry and is not 
simple to separate from transfer across thin walls, so both the 
rate and functional form of the transfer might differ from ex¬ 
periment. Prior work has shown that 2D CPM simulations of 
wet foams using the current method reproduce the bulk growth 
behaviors and distribution functions seen in experiments, in¬ 
cluding their dependence on foam wetness f\M . While such 
agreement shows that the functional form of the mass transfer 


via Plateau borders must be similar in CPM and wet foams, the 
agreement does not demonstrate that the functional forms are 
identical, nor how the rate depends on model parameters. In 
particular, for foams composed of a gas with fairly low solu¬ 
bility in the liquid phase, we expect gas transport via walls to 
be much faster than via Plateau borders. In our current CPM 
simulations, the interaction range is a significant fraction of the 
typical Plateau-border size, thus we expect the transport rates 
via walls and Plateau borders to be more similar and the differ¬ 
ence in growth rate between dry and wet foams to be less pro¬ 
nounced in the simulations than in an experiment with a foam 
with a low solubility gas phase. In the near future, we plan to 
examine the explicit relationship between Plateau border shape, 
size and model parameters, which will allow us to control the 
relative transport rates explicitly. 

We first relax the starting configuration for 200 MCS, to al¬ 
low the bubble and liquid volumes and geometries to adjust to 
generate a foam with the initially specified area distribution, 
then begin the simulation proper. To ensure that the equilibra¬ 
tion of wall and Plateau-border shapes is much faster than the 
rate of mass diffusion, after each 20 MCS of coarsening steps 
we run at least 10 MCS of relaxation steps, continuing until the 
number of accepted voxel copies during a relaxation MCS is 
zero or less than 1% of the number of accepted voxel copies 
during the previous relaxation MCS. Because mass transfer 
controls the time-scale of coarsening in experiment, with re¬ 
laxation occuring simultaneously, the simulation time t in our 
graphs reflects the number of coarsening steps only. During 
coarsening steps we set Agas = 0 so bubble volumes are free to 
vary and T = S for 0 = 0, and T = 25 for 0 > 0 since the 
volume constraint on the fluid particles increases their pinning 
to the grid. During relaxation steps, we set T = 0, Agas = 7, 
and for each bubble equal to the bubble’s volume at the 

beginning of the step. The volume constraint prevents mass 
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Figure 6: Bubble size distribution. Upper panel: Volume distributions as func¬ 
tion of V/ <y> for simulations starting from a log-normal bubble volume distri¬ 
bution and 0 = 0.20 (top). The black line shows the average over 10 simulation 
runs for 31 time points after the onset of self-similar growth. Lower panel: 
Self-similar volume distributions for different wetnesses. As for the growth 
rates, the distributions are similar but statistically distinct. 



Figure 7: Bubble face-number distribution. Upper panel: Face-number distri¬ 
butions for simulations starting from a log-normal bubble volume distribution 
for 0 = 0.20 (top). The black line shows the average over 10 simulation runs for 
31 time points after the onset of self-similar growth. Lower panel: Self-similar 
bubble face-number distributions for different wetnesses. 


3. Results 


transfer between bubbles, while allowing boundaries to move 
to minimize their interface energy. 

The simulations use the CompuCellSD environment ca, 
publicly available at http://www.compucell3d.org/ We ran the 
simulations on an Intel Core 17- 3700K processor. We ran 10 
replicas for each set of simulation parameters, starting from a 
log-normal volume distribution for each set of parameters, with 
0 = 0.00, 0.05, and 0.20. Each simulation replica took from 
2 to 7 days to run, depending on the number of simulated ob¬ 
jects (bubbles plus fluid particles), which varies with wetness. 
Because CompuCelBD limits the total number of objects in a 
simulation, we chose the size of the grid to respect that limit. 


Figure shows 2D sections of a simulated 3D foam be¬ 
fore (left) and after (right) a set of relaxation steps. Relax¬ 
ation smoothes the boundaries while preserving bubble vol¬ 
umes. Figure shows 3D images of a larger simulation for 
0 = 0.20, run on a 300^ grid: boundaries are smooth and the 
foam texture closely resembles that in the experiments in Figj^ 
To check whether the foam reached a self-similar growth 
regime we compared the bubble growth rate, G, as a function 
of the normalized bubble volume V/ {V) (flgurej^ at different 
times. The growth rates are identical within error for t > 600 
MCS for every wetness value (upper panel for 0 = 0.20, data 
not shown for 0 = 0.00,0.05). The lower panel shows the self¬ 
similar growth rates (averaged over 10 runs and 31 times) for 
0 = 0.00, 0.05, and 0.20. Figure S5 in the supplementary mate- 
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roughly one decade is compatible with a power law, which is 
a requirement for a self-similar growth regime |0. Fits for 
600 MCS < t < 3500 MCS yield volume scaling exponents of 
1.50, 1.33, and 1.20 + 0.01, for 0 = 0.00,0.05 and 0.20 respec¬ 
tively (corresponding to radius scaling exponents of 0.5, 0.44 
and 0.4). For t > 600 MCS, the number of bubble rearrange¬ 
ments Nti per MCS also decreases as a power of time (figure 
inset in the upper panel). The lower panel in figurej^ shows log- 
log plots of the decrease of the number of bubbles for 0 = 0.20, 
starting from initial configurations with normal, log-normal and 
bidisperse bubble volume distributions. The growth exponent is 
independent of the initial volume distribution. Figure S7 in the 
supplementary materials shows the self-similar distributions for 
number of faces for 0 = 0.20 and different initial distributions: 
the overlapping of these functions is a further evidence for the 


Figure 8: Bubble number vs time. Upper panel: Change in the number of 
bubbles for wetnesses (p = 0.00, 0.05 and 0.20, starting from a log-normal 
distribution of bubble volumes. After a transient period of roughly 600 MCS, 
the number of bubbles and the rate of Tls (inset) both decrease as a power of t. 
Lower panel: Dependence of the fractional decrease in the number of bubbles 
for different initial volume distributions with (p = 0.20. Averages and standard 
deviations calculated from 10 simulation replicas per data point. 


scaling regime. 


4. Discussion and conclusions 


rials shows a detail of this figure: the growth rates for 0 = 0.20 
significantly differ from those in the drier foam simulations. 

Figures and [ 7 ] show the distributions of normalized bubble 
volumes and face number. Both distributions are time invariant 
for t > 600 MCS for all wetnesses. The lower panels com¬ 
pare the average distributions for simulations with 0 = 0.00, 
0.05, and 0.20. The normalized volume distributions are simi¬ 
lar, though statistically distinct (see figure S6 in the supplemen¬ 
tary materials). The face-number distributions differ dramati¬ 
cally for foams with different wetnesses. 

The upper panel in figure shows log-log plots of the de¬ 
crease in number of bubbles starting from initial configurations 
with a log-normal distribution of bubble volumes. The straight 
line behavior observed for all wetnesses for t > 600 MCS over 


So far, experiments and simulations have investigated self¬ 
similar growth primarily in dry foams migiia. Together, the 
growth rates in figure and distribution functions in figures 
an d0 and the power laws seen in figure show that the foams 
reach a self-similar growth regime after a transient of 600 
MCS for all wetnesses and all three initial volume distributions. 

Radius scaling exponents in experimental foams depend on 
the wetness, dropping from 1/2 to 1/3 as 0 goes from zero 
to one. In our simulations, the radius scaling exponent for 
0 = 0.20 is significantly smaller than in the experiments of Isert 
ca but similar to that in Fortuna et al. ca. This difference 
between experiments, as well as between experiment and sim¬ 
ulation, may depend on the relative transport rates in thin walls 
and Plateau borders, which we will investigate in future work. 
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